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Statistical mechanics of the discrete nonlinear Schrodinger equation is studied by means of 
analytical and numerical techniques. The lower bound of the Hamiltonian permits the construction 
of standard Gibbsian equilibrium measures for positive temperatures. Beyond the line of T = oo, we 
identify a phase transition, through a discontinuity in the partition function. The phase transition 
is demonstrated to manifest itself in the creation of breather-like localized excitations. Interrelation 
between the statistical mechanics and the nonlinear dynamics of the system is explored numerically 
in both regimes. 
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The pioneering studies of Fermi, Pasta and Ulam ||] 
(FPU) showed that energy exchange between coupled 
systems may be suppressed in the presence of nonlin- 
earity; instead a type of behavior that severely contrasts 
equipartition among the hnear modes is observed. The 
question of whether equipartition of excitation energy al- 
ways appears is a contemporary issue in various fields 
of physics. Many manifestations of nonequilibrium and 
non-equipartition phenomena equivalent to the dynami- 
cal behavior of systems with few degrees of freedom con- 
trasting statistical mechanics expectations have been ob- 
served. Some of these phenomena, and therefore the ab- 
sence of immediate equipartition expressed in terms of 
self-trapping of energy, play an important role for opti- 
cal storage patterns in nonhncar fibers, condensed matter 
physics, and biophysics 

A particularity of discrete nonhnear systems is their 
ability to sustain strong localization of energy This 
is accomplished via intrinsic localized modes (breathers) 
which are modes that remain stable for extremely long 
times. So far it is a largely unaddressed problem how 
to handle and describe these excitations in a statistical 
mechanics framework although it has been argued that 
breathers may act as virtual bottlenecks Q delaying the 
thermalization process. 

In this work, we develop a statistical understanding of 
the dynamics, including the breathers, in a discrete non- 
linear Schrodinger (DNLS) equation. The DNLS equa- 
tion plays a significant role in several branches nonlinear 
physics as a simple physical model because it may ap- 
proximate many of the above mentioned nonlinear sys- 
tems. We study analytically and numerically the ther- 
malization of the lattice for T > 0. We identify the 
regime in phase space wherein regular statistical me- 
chanics considerations apply, and hence, thermalization 
is observed numerically and explored analytically using 
regular, grand-canonical, Gibbsian equilibrium measures. 



However, the nonlinear dynamics of the problem ren- 
ders permissible the realization of regimes of phase space 
which would formally correspond to "negative temper- 
atures" |] in the sense of statistical mechanics. The 
novel feature of these states is that the energy tends to 
be localized in certain lattice sites forming breather-like 
excitations. Returning to statistical mechanics, such re- 
alizations, which would formally correspond to negative 
temperatures, are not possible (since the Hamiltonian is 
unbounded from above, as is seen by a simple scaling 
argument similar to the continuum case ^) unless one 
refines the grand-canonical Gibbsian measure to correct 
for that. This correction will necessitate a discontinuity 
in the partition function signaling a phase transition that 
we identify, numerically, with the appearance of breather 
modes. 

In order to explore and illustrate the scenario described 
above we consider the one-dimensional DNLS equation in 
the form 
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(1) 



dH 



where the overdot denotes time derivative, m is a site in- 
dex, and V in a, tunable coefficient to the nonlinear term 
[Q . Equation (|^) is the equation of motion, 11)^ - 
where H. is the Hamiltonian function given by 



for which «V'mj V'm form canonically conjugate pairs of 
variables. In addition to the conserved energy 7i, the 
quantity A — IV'mPj is also conserved by the dy- 
namics of Eq. (^ and serves as the norm of the system. 

In order to study the statistical mechanics of the sys- 
tem, we calculate the classical grand-canonical partition 
function Z. We first apply the canonical transformation 
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Tpm = VA^ exp{i4>m) , leading to 



The partition function then becomes, 



Z= / l[dc^,ndA„,exp[~f3{n + f^M . (2) 

where the multipher fi is introduced in analogy with a 
chemical potential to ensure conservation of A. Straight- 
forward integration over the phase variable (j)m, reduces 
the symmetrized partition function to, 
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This integral can be evaluated exactly in the thermody- 
namic limit of a large system (iV —^ oo) using the eigen- 
functions and eigenvalues of the transfer integral operator 



dAm K{Am, Am+i) y{Am) = Xy{Am+i), 



where the kernel k is 



exp 



+ z^) + f^{x + z) 
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Similar calculations have been performed for the statisti- 
cal mechanics of the (j)^ field ^ , and for models of DNA 
denaturation One obtains Z ~ (27rAo)^, as N oo 
where Aq is the largest eigenvalue of the operator. From 
this expression the usual thermodynamic quantities such 
as the free energy, F, or specific heat can be calculated. 
More importantly, for our purpose we can obtain the av- 
eraged energy density, h = {Ti.)/N, and the average exci- 
tation norm, a ~ (A) /N as 
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The average excitation norm a can also be calculated 
as a = ;g- Hm ^^"i ^xp [-^ (H + ^A)] , where the 
integral again can be calculated using the transfer tech- 
nique and yields a = yQ{A)AdA, where j/o is 
the normalized eigenfunction corresponding to the largest 
eigenvalue Ao of the kernel k (Eq. (^)). This shows, that 
p{A) = VoiA) is the probability distribution function 
(PDF) for the amplitudes A. 

The problem is now reduced to finding the largest 
eigenvalue Aq and the corresponding eigenfunction i/q 



of the transfer operator, Eq. (^. This we do numeri- 
cally. However, two limits {P ^ oo and (3^0) are also 
amenable to analytical treatment. 

First, notice that the Hamiltonian is bounded from be- 
low and one can observe that this minimum is realized by 
a plane wave, Tpm — \/aeKpirmT, whose energy density is 
h = —2a + §a^. This relation defines zero temperature, 
or the (3 = oo line. 

The high temperature limit is also tractable. When 
/3 <C 1 the modified Bessel function in the transfer oper- 
ator can be approximated, to leading order, by unity (this 
amounts to neglecting the coupling term in the Hamilto- 
nian). This allows us to reduce the remaining eigenvalue 
problem to the approximate solution valid for thcrmal- 
ized independent units. 
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Using this approximation we can, enforcing the con- 
straint = 7 (where 7 remains finite as we take the 
limits /3 — > and /j, 00), obtain h = v/^"^ and a = I/7. 



Thus, we get h 



at /? = 0. 




FIG. 1. Parameter space {a,h), where the shaded area is 
inaccessible. The thick lines represents the /3 = 00 and (3 = 
lines and thus bound the Gibbsian regime. The dashed line 
represents the h = 2a + line along which the reported nu- 
merical simulations are performed (pointed by the symbols) . 

Figure 1, depicts (with thick lines) the two parabolas 
in (a, /i)-space corresponding to the T = and T = 00 
limits. Within this region all considerations of statistical 
mechanics in the grand-canonical ensemble are normally 
applicable and there is a one-to-one correspondence be- 
tween (a, /i) and (/?, /i). Thus, within this range of pa- 
rameter space one expects the system to thermalize in 
accordance with the Gibbsian formalism. However, the 
region of the parameter space that is experimentally (nu- 
merically) accessible is actually wider since it is possible 
to initialize the lattice at any energy density h and norm 
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density a above the T = line in an infinite system. 

A statistical treatment of the remaining domain of 
parameter space can be accomplished introducing for- 
mally negative temperatures. But the partition func- 
tion (^) is not suited for that purpose since the con- 
straint expressed in the grand-canonical form fails to 
bound the Hamiltonian from above. In all the alter- 
native approaches of the study of negative tempera- 
tures we will have to consider a finite system of size 
N. As suggested in wc can still consider the grand- 
canonical ensemble using the modified partition function 
Z'(/3,/x') = / eM-Pin + ^i'A^))l\,ndi'md^P*m, but this 
introduces long range coupling and ^' will have to be of 
order 1/N. Now /3 can be negative since H + fi'A^ can 
be seen to be bounded from above when fj,' < —v/2N. 
The important consequence of this explicit modification 
of the measure, is a jump discontinuity in the partition 
function, that in turn signals a phase transition. More 
explicitly, if one starts in a positive T, thermalizable (in 
the Gibbsian sense) state in phase space with /i > 0, and 
continuously varies the norm, then one will, inevitably, 
encounter the /3 = parabola. Hence, in order to proceed 
in a continuous way, a discontinuity has to be assigned 
to the chemical potential. This discontinuity will destroy 
the analyticity of the partition function as the transition 
line is crossed, and will indicate a phase transformation 
according to standard statistical mechanics. 

From the microcanonical point of view it is also natu- 
ral to consider negative temperatures because it is possi- 
ble to maximize the energy under the constraint of fixed 
norm in a finite system. It can be seen that the configu- 
ration which realizes this maximum is an exact breather 
solution, whose total energy and frequency scale as 
and N , respectively. Thus, the number of microstates 
sharing the same energy E will decrease with increasing 
E if the norm A is kept fixed. Due to the definition of 
temperature (1/T — dS/dE\^, T becomes negative at 
high energy density and the [3 = Q line is the line of max- 
imum entropy. Actually we can say that the constraint 
of fixed norm ^ is a "topological" reason for large ampli- 
tude breather-like excitations to be expected to appear. 

In order to characterize the dynamics of both phases 
(above and under the /3 = line) and to verify that 
the system does relax to a thermalized state, we per- 
form numerical experiments. We restrict the parameters 
(a, h) to the dashed-line of Fig. |l|, choosing a perturbed 
phonon with wavevector q = Q {ipm = v^)i which 
the energy norm relationship is ft. = 2a -I- as initial 
condition. An infinitesimal perturbation to such a lin- 
early unstable mode for systems of oscillators or nonlin- 
ear ly coupled particles, is well-known [^|o| to give rise 
to long-lived localized excitations via modulational in- 
stabilities. For these initial conditions, the important 
question is whether the same phenomenology appears in 
the DNLS system; i.e., whether relaxation to equilibrium 
is really achieved and whether we can observe different 



qualitative behavior on the two sides of the /3 = line. 

Figure 2 shows three typical examples of what can be 
observed when the energy-norm density point lies below 
the /3 = line (the symbols refer to Fig. 1). The q — 
wave is unstable and the energy density forms small lo- 
calized excitations but their lifetime is not very long and, 
rapidly, a stationary distribution of the amplitudes 
is reached (Fig. 2). Different kinds of initial conditions 
(with same energy and norm densities) produced the 
same results. In conclusion, the system reaches an equi- 
librium state which is perfectly recovered by means of the 
transfer operator method. Moreover it can be checked on 
Fig. 2 that the curvature of \ogp{A) (i.e., —(i) tends to 
zero when h = va?' . (The cut-off at high amplitudes is 
due to finite size effects). In this domain of parameter 
space, high amplitude excitations are highly improbable 
and can be considered as simple fluctuations; as shown on 
Fig. ^, large amplitude fluctuations have been recorded 
but were checked in the evolution pattern to disappear 
rapidly. 




FIG. 2. Distribution oi A — \tp'\^ for three cases under (and 
on) tlie transition line. The solid lines show the results of sim- 
ulations and the symbols are given by the transfer operator. 
Curves are vertically shifted to facilitate visualization. 



The scenario is very different when the energy and 
norm densities are above the /3 = line. We can ob- 
serve a rapid creation of breather excitations due to the 
modulational instability accompanied by thermalization 
of the rest of the lattice. Once created, these localized ex- 
citations remain mostly pinned and because the internal 
frequency increases with amplitude their coupling with 
the small amplitude radiation is very small. This intro- 
duces a new time scale in the thermalization process ne- 
cessitating simplectic integration for as long as 10^ — 10'' 
time units in order to reach a stationary PDF. This can 
also be qualitatively justified by the effective long range 
interactions, introduced in the modified partition func- 
tion, which will produce stronger memory effects as one 
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observes regimes in phase space which are further away 
from the transition hne (since the long range interaction 
will be stronger). 

Typical distribution functions of the amplitudes are 
shown in Fig. ^. The presence of high amplitude exci- 
tations is directly seen here (more straightforwardly we 
observed standing breathers in the spatial pattern) . The 
dotted line represents the PDF in the case where the ini- 
tial condition is chosen at random, using a larger system 
size: we check that the initial condition seems unimpor- 
tant, but the system size does influence the amplitude 
of the highest breathers. The cut-off value in the very 
large system limit, as well as the persistence of a bump 
in the PDF, is still an open question, since we have no 
prediction for the PDF above the (3 = line. However the 
positive curvature of the PDF at small amplitudes clearly 
indicates that the system evolves in a regime of negative 
temperature and the appearance of the phase transfor- 
mation is signified in the dynamics by the appearance of 
these strongly localized persistent breathers. 



ified measure to ensure normalizability, which therefore 
necessitates a phase transition leading to these localized 
modes. Our numerical simulations strongly support this 
theoretical picture, illuminating this novel quantitative 
connection between nonlinear dynamics and statistical 
mechanics. 
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FIG. 3. Distribution oi A = \^p\ for parameters (h,a) 
above the transition line (triangles and stars as in Fig. 
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The actual dynamics in the negative temperature 
regime is studied more closely in Ref. pT| . 

Finally, we can draw an interesting parallel with what 
has been known in plasma physics and hydrodynamics 
for several years , where the appearance of localized 
structures (of vortices in that case) is also related to a 
description in terms of negative temperatures. 

In conclusion, studying the DNLS system, we have 
been able to quantify and explain, through analytical 
calculations supported by numerical computations, the 
behavior in different regimes of the [h, a) phase space. 
We have been able to link the regime of thermalization 
to the regime where regular statistical mechanics is ap- 
plicable in the Gibbsian sense. Further, we have traced 
the explanation of the appearance of localized modes in 
different regimes in phase space to the need for a mod- 
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